Simulability and regularity of complex quantum systems 
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We show that the transition from regular to chaotic spectral statistics in interacting many-body 
quantum systems has an unambiguous signature in the distribution of Schmidt coefficients dynami- 
cally generated from a generic initial state, and thus limits the efficiency of the t-DMRG algorithm. 
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Complexity is paradigmatic in many areas not alone 
in physics, but equally so in the life and social sciences 
and in economics The characteristic property of a 
"complex system" resides in the difficulty of its efficient 
simulation through reduction to manageable size. More 
formally, e.g. the minimum length of an algorithm de- 
signed to simulate the system under study can serve as a 
quantitative measure of complexity ■ Complexity thus 
has its very tangible counterpart in the numerical over- 
head required for an accurate simulation, and implies a 
challenge for computational physics: any time we suc- 
ceed to minimize that overhead we actually prove that 
the underlying complexity is smaller than anticipated. 

Complex systems abound in nature, from interacting 
many-particle systems to deterministic chaos in few de- 
grees of freedom 0] , and from macroscopic to microscopic 
scales. On the quantum level, complexity has its counter- 
part in complex spectral structures, described by random 
matrix theory and by now well-established in nuclear [J] 
and atomic physics @, mesoscopics @, and photonics 
0, H[ • Recently, "designed complexity" moved into reach 
for state of the art experiments on ultracold atoms in 
periodic optical potentials, and considerable effort is de- 
voted to implementing solid state Hamiltonians with un- 
precedented control in such systems 0, [l(| EL E3 ■ Since 
analytical treatments are often unavailable to describe 
such many-particle dynamics, efficient numerical tools 
are in need, and "efficient" means here that the required 
numerical resources as memory and/or execution time 
scale favorably as compared to the exponential growth of 
Hilbcrt space dimension with system size. 

For ID systems, the advent of Density Matrix Renor- 
malization Group methods [H, has aided efficient 
simulation considerably. There the underlying idea lies 
in the construction of a suitable local basis that makes 
it possible to represent a system state in terms of sig- 
nificantly fewer basis states than the total dimension of 
Hilbert space suggests. More specifically, this approach 
makes use of a truncated Matrix Product State (MPS) 
ansatz 15], to reduce the number of coefficients required 
to specify the state to a manageable number. In partic- 
ular in perturbative regimes, where a system has a nat- 
ural basis, such techniques work very successfully, and 
the ground states of ID systems with local Hamiltoni- 



ans are typically well represented in this form [l6|, [T3|- 
Also dynamics can be tackled by time-dependent Density 
Matrix Renormalization Group techniques (t-DMRG). 
These methods have been shown to work well for low- 
energy initial states [H, EM 13, giving rise to near- 
equilibrium time-evolution. For a generic initial state, 
however, they may work only for short times [20(, and ap- 
parently depend strongly on the properties of the Hamil- 
tonian. 

Here we investigate the connection between complexity 
and t-DMRG methods: Can t-DMRG efficiently simulate 
complex many-particle systems in general, or does the 
efficiency of t-DMRG rather identify parameter regimes 
where the dynamics - governed by the underlying spec- 
tral structure - is actually rather "regular"? To assess 
this question, we need independent measures of complex- 
ity, which, on the spectral level, are provided by the the- 
ory of quantum chaos Q. As we show in the following, 
the appearance of "globally irregular" vulgo "chaotic" 
spectral structure implies the breakdown of simulability 
by t-DMRG techniques. More specifically, we observe a 
characteristic qualitative change of the distribution of dy- 
namically generated Schmidt coefficients that directly re- 
flects the defining property of complexity, namely the ef- 
fective irreducibility of the Hilbert space dimension: Any 
basis truncation leads to the rapid accumulation of un- 
controllable errors in the simulation. 

Specifically, we investigate the tilted Bose-Hubbard 
model. Since our results rely on universal statistical 
properties that do not depend on system specificities, the 
observed phenomena are expected to hold more generally. 

The tilted Bose-Hubbard model we consider here cor- 
responds to bosonic atoms trapped in the lowest band 
of an optical lattice subject to an additional static field. 
For m lattice sites, the system is described by the Hamil- 
tonian 
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where a] (a;) are creation (annihilation) operators for a 

particle at lattice site I, and hi = a\ai is the correspond- 
ing number operator. The system's dynamical properties 
are characterized by the tunneling constant J, the inter- 
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action strength U, and a linear potential with strength 
F, resulting, e.g., from gravity when the optical lattice 
is tilted with respect to gravitational equipotential lines. 
The rich dynamics of this system ranges from perfect os- 
cillations [2ll[2l| with the Bloch period T B = 2n/F, for 
weak interactions, to chaotic dynamics for situations in 
which interaction and tunneling are of comparable order 
of magnitude [231 ]. 

As mentioned above, simulability via t-DMRG meth- 
ods relies on the possibility to effectively decimate the 
Hilbert space of the system in the course of the propa- 
gation, by using a truncated MPS representation. This 
is based on the Schmidt decomposition for every possible 
bipartite splitting of the system [2^ |. 

IV>HEa q |$W)|$L b ]> . (2) 

a=l 

where A and B label the two subsystems with Schmidt 
eigenstates {l^k 4 ')} and {|$^}}, and the real coefficients 
A Q , with normalization condition J2a^a = 1; describe 
the quantum correlations between the two subsystems. 
In the truncated MPS, we make an approximation to the 
state by setting an upper bound \ on Xab , retaining only 
those eigenstates with the largest X a . This approxima- 
tion is good if the X a decay rapidly as a function of the 
index a, when ordered from largest to smallest. Thus, it 
relies on the assumption that the entanglement between 
two parts of the system, as, e.g., measured by the von 
Neumann entropy 

S=-J2^og 2 X 2 a , (3) 

a 

is never too large. The von Neumann entropy thus also 
provides an estimate for the value of x required to rep- 
resent the state, and its behaviour is often used as an in- 
dicator of simulability of the system. If S grows rapidly 
as a function of time, simulation of the system will be 
difficult, whereas if S is bounded during the dynamics, 
then we can fix x and compute the dynamics over long 
time periods. 

We now consider this approximation for the time evo- 
lution of various initial states, in different parameter 
regimes of the Hamiltonian, Eq. |T]), above. In the fol- 
lowing we always discuss that bipartition that yields the 
largest von Neumann entropy, since this limits the effi- 
ciency of the t-DMRG algorithm (However, other bipar- 
titions also yield qualitatively similar results.). In order 
to permit our subsequent comparison with the system's 
spectral properties, we choose a relatively small system, 
beginning with eight bosons in eight neighboring sites at 
the centre of a lattice of length 64, with Dirichlet bound- 
ary conditions, for different values of U / J and F/J. Note 
that for F I J > 1, the particles hardly spread on the time 
scale of the simulation, which makes it possible to com- 
pare the dynamical behavior of the larger system to the 
spectrum of the Floquet-Bloch operator on nine sites. 
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FIG. 1: Distribution of average Schmidt coefficients \ a , ob- 
tained from averaging over ten initial states of eight particles 
initially distributed over eight lattice sites. The total grid size 
was 64, to eliminate boundary effects. Shown are the x — 100 
largest average Schmidt coefficients, sorted in descending or- 
der, for U/J = 1 (thick lines) and U/J = 10 (thin lines), and 
for different values of the static tilt F/J = 1.0 (solid), 1.5 
(dashed), 2.0 (dotted), at t = 4.776 x T B . 



However, for very weak static fields the atoms rather dif- 
fuse through the system, and travel into initially unoccu- 
pied regions, so that the correspondence between spectral 
statistics and t-DMRG is lost. Therefore, our subsequent 
analysis does not extend to very small values of F/J, 
where diffusion sets in. 

Fig. [T] shows the averaged Schmidt coefficients af- 
ter an evolution time t = 4.776 x Tg, with the aver- 
age taken over ten separable initial states of the form 
iV'init) ~ \ni) <g> \n,2) <8> • ■ • <8> \n%) (each with a different 
random realisation of occupation numbers ni of eight ini- 
tially occupied sites). The distributions can be divided 
into two qualitatively different categories: For strong in- 
teraction U/J = 10, few comparably large coefficients 
dominate the distribution which exhibits a rapidly di- 
minishing tail. In contrast, for weak interaction U / J = 1, 
the distribution shows a slowly decaying tail with many 
non-negligible coefficients of comparable weight. 

The eligibility of two such distributions for the basis 
truncation required to apply the t-DMRG algorithm is 
fundamentally different: while a distribution of the for- 
mer type allows the truncation of the major portion of 
the Schmidt basis with virtually vanishing loss of accu- 
racy, dropping a few basis states in the latter case already 
will lead to a sizable error. That is, we can identify situa- 
tions that can be described efficiently with t-DMRG, but 
there are parameter regimes of the same system where a 
faithful representation of the solution must be spanned 
essentially by the complete Hilbert space. Any numerical 
simulation then is plagued by highly unfavorable scal- 
ing. In particular, as immediately evident from Fig. [TJ 
MPS basis truncation at x = 100 for U/J = 1, F/J = 1 
enforces a rapid decrease of the Schmidt coefficients for 
a > 80, which, in turn, will induce artifacts in the sim- 
ulated dynamics (as we confirmed by running computa- 
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tions for varying x). 

Such transition from efficient to inefficient representa- 
tion in a MPS basis has its cause in a sudden and pro- 
nounced transition in the underlying spectral structure, 
as we will now evidence by direct inspection of the spec- 
trum of the time evolution (or Floquet) operator 

U(T B ) = fexp ^~iJ o H(i)dtj (4) 

generated by the time-dependent, transformed Hamilto- 
nian 

1=1 i=i 

(T denotes time ordering). Due to the translational 
invariance of H(t) with periodic boundary conditions, 
U(Tb) decomposes into the direct sum of operators la- 
beled by distinct values of quasimomentum n (23[, and 
so does its spectrum. The statistical analysis therefore 
requires a diagonalization at fixed quasimomentum, and 
we chose n = here. 

The integrated level spacing distribution I(s) of the 
eigen- frequencies of U(Tb) allows to distinguish regular 
spectral structure (tantamount of weakly coupled basis 
states) , described by Poissonian statistics [25| 

I P (s) = [ P(s')ds' = 1 - exp (-s) , (6) 
Jo 

from a chaotic spectrum that obeys Wigner-Dyson statis- 
tics [4j, 

I w {s) = f-cxp(-^ 2 /4) . (7) 
Fig. [5] displays the mean square deviation 

A2 = JfJ ds(f(s)-I(s)) 2 (8) 

of the numerically obtained spectra with respect to Iw(s) 
and Ip(s). 

For U/J — 10, the system obeys Poissonian statistics, 
irrespective of F/J. However, in the case of U/J = 1 
there are three different regimes: for F/J < f.3, the de- 
viation from (irregular/chaotic) Wigner-Dyson statistics 
is negligible, for f.3 < F/J < 2 the distribution changes 
its character, and turns (regular) Poissonian for F/J > 2. 
This transition is also reflected in the entropy of Schmidt 
coefficients S, Eq. ([3]), that is depicted as a function of 
F/J after a simulation time of t = 4.776 x Tp in Fig. [3] 
The inset shows the entropy as a function of time, for 
different values of the static field strength F/J. In the 
regime of regular level statistics, S grows only initially, 
whereas it keeps increasing in the chaotic regime. As 
a matter of fact, the saturation of S for F/J < 1 and 
U / J — 1 is once again a numerical artifact of the trunca- 
tion at x = 100 ; The simulation only approximates the 




FIG. 2: Mean square deviation A 2 , Eq. ([SJl, of the Floquet- 
Bloch operator's Q nearest neighbor distribution from Pois- 
sonian (dashed lines) and Wigner-Dyson statistics (solid 
lines), for eight particles on nine lattice sites, as a function of 
F/J. U/J = 1 (thick lines), 10 (thin lines). 
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FIG. 3: Many-particle entanglement vs. F/J, measured by 
the average von Neumann entropy S, extracted from t-DMRG 
simulations of eight particles for ten different initial states at 
t = 4.776 x T B . U/J = 1 (thick), 10 (thin line). Inset: S vs. 
t, for various tilt strengths F/J = 1.0 (solid), 1.5 (dashed), 
2.0 (dotted lines), and the same values of U/J. 
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FIG. 4: Average number of Schmidt coefficients larger than 
e = 0.01, vs. F/J, after a simulation time t — 4.776 x Tb- 
U/J = 1 (thick), 10 (thin lines). 
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real values of S from below, i.e., at given t, S grows when 
increasing \. This behavior of S is perfectly consistent 
with our observations on Fig. [1] Note that the slight in- 
crease of S towards small values of F/J, for U/J = 10 
in Fig. [3j corresponds to a narrow distribution of the 
Schmidt coefficients in Fig.[IJ and therefore does not hin- 
der efficient simulation - in perfect agreement with the 
regular spectral structure in this regime spelled out by 
Fig.H 

While the signature of the sharp "chaos-transition" ob- 
served in Fig.[5]is somewhat smoother in the correspond- 
ing behaviour of the von Neumann entropy in Fig. [3J 
the number of Schmidt coefficients larger than a cer- 
tain threshold e turns out to provide an equally sensitive 
probe as the spectral statistics, as demonstrated in Fig. [4] 
(for e = 0.01): Whereas in the regular regime (U/J= 10, 
or U/J = 1 with F/J > 2) less than 20% of the co- 
efficients exceed the threshold e, essentially all of them 
contribute in the chaotic regime (U/J = 1, F/J < 1.3). 
That is, whereas the t-DMRG algorithm allows an effi- 
cient simulation of the Bose-Hubbard dynamics in the 
regular regime, a basis truncation in the chaotic regime 
will rapidly lead to sizable errors in the simulation. An 
accurate description requires large numerical efforts that 
scale exponentially, much as the system size itself. This 



observation also holds for larger systems, where t-DMRG 
is a powerful tool in the regime of regular spectral struc- 
ture, and where an exact treatment of the dynamics be- 
comes unfeasible: In t-DMRG calculations with 20 atoms 
in 20 sites we observe precisely the same characteristic 
changes in the distribution of Schmidt coefficients as ob- 
served in Figs. [T] and 21 Therefore, the inefficiency of 
t-DMRG simulations, quantified by the statistical quan- 
tities evaluated in Figs. [2] and [4] is an unambiguous indi- 
cator of the underlying complexity of the many-particle 
dynamics. 

Moreover, since the non-existence of a natural basis is 
a rather generic feature of quantum chaotic systems, we 
expect that our observations directly translate to generic 
many-body quantum systems. In fact, we conjecture that 
the distributions of Schmidt coefficients of typical states 
in spectrally re gula r and chaotic systems also exhibit uni- 
versal features [26| - much like the energy level distribu- 
tions of regular and chaotic quantum systems. 
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